Structure prediction based on ab initio simulated annealing for boron nitride 
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Possible crystalline modifications of chemical compounds at low temperatures correspond to local 
minima of the energy landscape. Determining these minima via simulated annealing is one method 
for the prediction of crystal structures, where the number of atoms per unit cell is the only informa- 
tion used. It is demonstrated that this method can be applied to covalent systems, at the example 
of boron nitride, using ab initio energies in all stages of the optimization, i.e. both during the global 
search and the subsequent local optimization. Ten low lying structure candidates are presented, 
including both layered structures and 3d-network structures such as the wurtzite and zincblende 
types, as well as a structure corresponding to the /3-BeO type. 

PACS numbers: 61.50.Ah 71.15.Nc 61.66.Fn 



I. INTRODUCTION 

The knowledge of the crystal structure of a solid 
compound is one of the basic questions in solid state 
theory [J, S S II H- Since the early 1990's, an effort 
has been made to develop methods to predict structures 
of solids without any experimental information about the 
structure. The starting point is the realization that any 
(meta)stable modification of a (solid) compound corre- 
sponds to a locally ergodic region on the energy land- 
scape of the chemical system. For low temperatures, such 
regions are centered on local minima of the energy that 
are surrounded by sufficiently high energy barriers. Since 
both the thermodynamically stable and the multitude of 
kinetically stable modifications are of interest, the global 
search is not restricted to the determination of the global 
minimum, but also includes local minima [J]. 

The most common methods used for the structure pre- 
diction of solids are simulated annealing [a LP , JaK genetic 
algorithms [1, Qil [HI, HH , basin hopping [lj [l4| , or the 
recently introduced metadynamics [la ]. 

Structure prediction usually involves a huge amount of 
CPU time, and therefore efficient ways to keep the calcu- 
lations tractable have to be found. Thus, the procedures 
were initially split in two steps: first, a global search on 
the potential surface was performed. The energy was 
evaluated with empirical potentials, e.g. Coulomb and 
Lennard- Jones potentials, or chemically/physically mo- 
tivated cost functions. After the global search, e.g. using 
simulated annealing, possible structure candidates were 
locally optimized on the ab initio level, usually with den- 
sity functional theory. 

Empirical potentials are very efficient, but also have 
various drawbacks: they work reasonably well for ionic 
systems, but less for covalent systems; and some knowl- 
edge of the expected bond type is required in advance 
to choose the potential. Recently [la] , we demonstrated 
that a full ab initio treatment is feasible in both stages, 
i.e. the global search and the subsequent local optimiza- 
tion can both be performed on the ab initio level. The 
system considered (lithium fluoride) was chosen for sev- 
eral reasons: the small number of electrons leads to fast 



calculations, and the ionicity of the system makes con- 
vergence easy. Moreover, the system had earlier been 
studied with model potentials [171 UM > an d it turned out 
that the relevant minima were the same when full ab ini- 
tio structure prediction was performed [l6| (for a brief 
summary see also (l9| ) . 

In the present work, boron nitride (BN) was chosen as 
an example for a covalent system. This is a significant 
extension of the previous work, since a covalent system 
is much more difficult to study: covalent bonds between 
the atoms have to be formed, and convergence for ran- 
dom structures is much more difficult in this case. At 
zero pressure, BN has a hexagonal structure, with space 
group 194, see e.g. 20]. Under pressure, it may trans- 
form to a zincblende 21] or a wurtzite structure 22j, and 
a corresponding phase diagram was obtained [23 1 . There 
is, however, an ongoing discussion concerning the cor- 
rect phase diagram, see e.g. [24} . Also, a rhombohedral 
structure was found (a layered structure similar to the 
hexagonal structure) [251 ] , and turbostratic boron nitride 
with a random arrangement of the layers has been re- 
ported [26j] . These structures are summarized in table [J 
For overviews, see also [271 128I]. Early ab-initio calcula- 
tions were performed from the mid 1980's onwards, e.g. 

M, M [MM HI, M M HE S3. 



The present task is for several reasons far from 
straightforward: first, the CPU time has to be reduced 
by a large factor. As we showed for the system LiF [Taj , 
simply performing simulated annealing with standard ab 
initio calculations would lead to CPU times of the order 
of 2 years for a single run, and often hundreds of such 
runs need to be performed. Therefore, very subtle meth- 
ods have to be employed to reduce the CPU time for the 
ab initio calculations, without losing too much accuracy. 
In our earlier study of LiF, very careful tests were re- 
quired in order to reduce the CPU time to a few days. 
Similarly, many tests have to be performed to reduce the 
required CPU time in the case of BN to reasonable values. 
Secondly, one needs a strategy to converge the system in 
all steps of the global exploration: at the beginning of 
the search, the unit cell has a very large volume, and the 
atoms are at random positions within the cell, i.e. the 
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TABLE I: The experimentally found structures. 



structure type cell parameters and fractional coordinates 
hexagonal BN [20] a=2.504 A, c=6.661 A 
B (1/3, 2/3, 1/4) 
N (1/3, 2/3, 3/4) 

zincblende [211 a=3.615 A 

B (0, 0, 0) 
N (1/4, 1/4, 1/4) 

wurtzitc [22] a=2.55 A, c=4.20 A 

rhombohedral [25] a=2.52 A, c=10.02 A 



configuration is like in a gaseous state. The total en- 
ergy calculation at such a geometry has to be converged, 
and it has to be converged for all the subsequent steps of 
the simulated annealing procedure. This has to be done 
in an automatic way, since thousands of calculations at 
successive geometries are performed. Finally, it is a very 
interesting and important question whether this proce- 
dure will find all the very different structure types: the 
layered structures, and the three-dimensional structures 
such as wurtzitc and zincblende. 



II. METHOD 

The general method consists of several steps: first, 
a simulated annealing run with a subsequent stochastic 
quench is performed, to identify possible candidate struc- 
tures. This is followed by a local optimization based on 
analytical gradients. Finally, the symmetry and space 
group are identified. This is repeated many times, in or- 
der to identify as large as possible a set of structure candi- 
dates, and to obtain some statistics about the structures 
found. 

The details for the BN calculations are described in 
the following paragraphs. 4 boron and 4 nitrogen atoms 
were placed at random positions in a large unit cell which 
was initially cubic with the cell parameter a=6.23 A. 
This initial volume is computed by first employing the 
atomic/ionic radii to estimate the total volume occupied 
by the atoms in the solid state and then multiplying this 
number by a sufficiently large factor (typically between 
3 and 10), to allow the atoms enough freedom to reach 
any arrangement independent of the initial random place- 
ment in the simulation cell. This factor was varied in 
preliminary calculations, where it was found that choos- 
ing a larger volume than the one selected (6.23 3 A 3 = 
242 A 3 , which is already by a factor of ~5 larger than 
the volume of e.g. the zincblende structure) did not lead 



to significantly different results. 

Each simulated annealing run had a length of 12500 
steps, and the initial temperature of 1.00 eV (correspond- 
ing to 11604 Kelvin) was reduced to ~ 0.78 eV at the end 
of the run. The length is thus somewhat longer than in 
the case of LiF [l6j], because it is more difficult to ap- 
proach possible candidate structures in a system with 
covalent bonds. The simulated annealing was followed 
by a quench with 5000 steps, i.e. a simulated anneal- 
ing run with a temperature of eV, which means that 
only downhill moves are allowed during the quench. The 
quench thus moves the geometry obtained after the sim- 
ulated annealing further towards a local minimum. The 
moves were chosen as: moving individual atoms (70%), 
exchanging atoms (10%), changing the lattice parame- 
ters with fixed fractional coordinates (10%), changing the 
lattice parameters with fixed cartesian coordinates (5%), 
and changing the origin (5%, this move is important if 
subsequently the cell parameters change). No symme- 
try was prescribed during the simulated annealing and 
quench runs, i.e. the space group was always PI. 

A minimum distance between two atoms (given by the 
sum of the radii of the atoms, multiplied with 0.7) was 
prescribed in order to avoid unrealistic geometries which 
may lead to numerical instabilities. The radii used were 
based on tabulated values for atomic and ionic radii, as a 
function of charge, and the Mulliken charge computed for 
the previous configuration. In those moves which change 
the lattice constant, the probability of reducing the lat- 
tice constant was enlarged to 60%, to speed up the re- 
duction of the cell size. 

The ab initio calculations were performed with the 
CRYSTAL06 code 3$] . which is based on local Gaus- 
sian type orbitals. Two basis sets were used during the 
simulated annealing runs, starting from a [3s2p] basis set 
for B and N, with the inner [2slp] exponents as in [39] ]. 
In the case of basis set I, additional sp exponents of 0.4 
for B and 0.3 for N were chosen, and the outermost ex- 
ponent of the B 2sp contraction (0.4652) was removed. 
In the case of basis set II, the outer sp exponents were 
chosen as 0.25 for B and 0.297 for N. The basis sets are 
given in table IIIIl In the stage of the local optimization, 
the basis sets used were the [3s2pl<i] basis sets from [32[ 
(basis set III in table IIII[) . 

The basis sets during the simulated annealing are 
therefore chosen slightly different from the ones used in 
the local optimization (less diffuse functions, no polar- 
ization functions), to speed up the global search which 
is the time-consuming part of the procedure. To test 
these basis sets, as a preliminary step, the energies of 
the wurtzite and of the layered Bk structure were com- 
puted with various basis sets (table [TTJ) . With basis set 
III, which is used during the local optimization, the Bk 
structure is more favorable by ~ 30 mEh, i.e. ~ 0.8 eV 
(lE h = 1 hartree = 27.2114 eV). The smallest basis set I 
gives preference to the wurtzite structure instead, by 33 
mEh (0.9 eV), and basis set II gives preference to the Bk 
structure, by 8 mEh (0.2 eV). One might thus fear that 
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TABLE II: The Hartree-Fock energies (in hartree, per four 
formula units) of the Bk and wurtzite structures, computed 
with the small basis sets used during the global optimization 
(basis set I, II), in comparison with the energy obtained with 
the basis set used for the local optimization (basis set III). The 
geometry was fixed at the computed equilibrium geometry of 
basis set III (table N} . 



structure basis set total energy 

B k I -316.6687 

II -316.7495 

III -316.8753 
wurtzite I -316.7014 

II -316.7419 

III -316.8452 



basis set I would not yield the layered structures during 
the global search; however, this is not the will 
be shown in the results section; and there does not seem 
to be a strong bias due to the basis set. Basis set I is 
however advantageous because calculations are roughly 
twice as fast as with basis set II. 

Concerning the choice of the ab initio method, it has to 
be taken into account that convergence at a random ge- 
ometry is necessary during the global exploration stage. 
As mentioned earlier, the initial geometry has a cell vol- 
ume ~ 5 times larger than the experimental volume, and 
the atoms are randomly arranged. A typical band struc- 
ture of such a geometry is very localized due to the large 
interatomic distances, and completely different from the 
band structure of the experimental geometry in the solid 
state. Therefore, convergence is absolutely non-trivial. It 
turned out that convergence was best achieved with the 
Hartree-Fock approach, due to the fact that the band 
gaps are usually very large. Indeed, the band structure 
and corresponding densities of states display band gaps of 
~ 6 eV (Hartree-Fock), ~ 0.5 eV (B3LYP), and - 0.1 eV 
(LDA) for this initial structure. Note that this gap corre- 
sponds to a random initial structure and is very different 
from the gap of the final structure; but it is necessary to 
converge a calculation for this random initial geometry, 
and for all the geometries subsequently generated, until 
the end of the simulated annealing and quench. For com- 
parison, calculations with the hybrid functional B3LYP 
were found to be much more difficult to converge, and 
a large mixing ratio was required: 90%, in combination 
with the Anderson mixing scheme; 35% was sufficient in 
the case of Hartree-Fock (the mixing ratio is the ratio 
with which the previous Fock operator is added to the 
new one, in order to achieve convergence). This leads to 
many iterations and thus a large CPU time. The local 
density approximation (LDA) was very difficult to con- 
verge for random atom arrangements, and needed more 
k points and level shifting. 



Interestingly, even for the initial structure consisting of 
widely spaced, nearly isolated atoms, it was sufficient to 
use the restricted Hartree-Fock approach, i.e. it was not 
necessary to take into account spin-polarization which 
would be the case for a free atom. 

The thresholds for integral selection were enlarged 
from 10~ 6 , 10~ 6 , 10~ 6 , 10" 6 , 10~ 12 to 10~ 4 , 10~ 4 , 10~ 4 , 
10~ 4 , 10~ 8 , respectively, and the self consistent field cy- 
cles were stopped when the difference between two sub- 
sequent cycles was below 10 -4 Eh- A mesh with 4x4x4 
fc-points was used. The error associated with the k mesh 
can be estimated by computing the energy difference 
when changing the lattice constant: e.g. changing the 
lattice constant of BN in the zincblende structure from 
3.7 to 3.6 A changes the energy by 0.01368 Eh for four for- 
mula units with a 4x4x4 mesh, and by 0.01436 Eh with 
a 8x8x8 mesh. The associated error with the mesh is 
thus 0.01436-0.01368=0.00068 E h and reasonably small. 

The simulated annealing and subsequent quench arc 
the time consuming parts, and a single run takes of the 
order of one week on a single CPU. The same approach 
would have been feasible with the B3LYP functional, but 
at a much higher cost, for the reasons discussed above 
(around a month instead of one week CPU time). It ap- 
peared therefore more reasonable to perform four times 
as many runs (here: around 329, see table HVj) using 
the Hartree-Fock approach, as with the B3LYP approach 
where around 80 runs would have been feasible with a 
comparable total CPU time. 

The local optimization is not very time consuming and 
was done with default parameters for the integral selec- 
tion and the self consistent field cycles. The full geometry 
optimization can by now be routinely performed with an- 
alytical gradients O E2, IS IS EH as implemented in 
the CRYSTAL06 release. The local optimization, start- 
ing from the structure after the quench, was done both at 
the HF level and at the LDA level; in nearly all the cases, 
the resulting final minimum structures turned out to be 
the same. In addition, for these final structures, also a 
B3LYP optimization was performed, in order to compare 
Hartree-Fock, B3LYP and LDA. The basis set used was 
basis set III in tabic IIII1 In addition, in appendix [XJ 
larger basis sets were tested, for comparison. 

The symmetry was analyzed with the program 
KPLOT [45| where algorithms to find the symmetry and 
space group [S |4?| are implemented. 

For the most important structures, the enthalpy was 
computed, in order to investigate the pressure depen- 
dence of the phases. 



III. RESULTS AND DISCUSSION 

The most relevant structures found are dis play ed in 
figures 1112131 and [4] visualized with XCrysDen [48| . Op- 
timized geometries are given in table [Vl First, the exper- 
imentally observed Bk structure was obtained (so-called 
hexagonal boron nitride, space group 194). Closely re- 
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lated are two additional layered structures, with space 
group 160 and 187, respectively. In the B k structure, 
sheets are made of edge-connected six-membered rings 
of 3 boron and 3 nitrogen atoms in alternating sequence 
(see figure [1] left). The neighboring sheets are stacked 
vertically below and above, with alternating atoms (i.e. 
N sits vertically above B, and vice versa; the stacking 
order is ABAB). In space group 160, the same sheets are 
formed, but only three atoms have neighbors in the layer 
above, and the other three atoms have neighbors in the 
layer below (the stacking order is ABCABC, i.e. like in 
rhombohedral BN [HI). The structure with space group 
187 has stacking order ABAB, where again three atoms 
have neighbors vertically above and below in the next 
layer. These three layered structures have a very simi- 
lar total energy, and also the enthalpy as a function of 
pressure looks very similar (see figure [5]). 

The wurtzite structure is displayed in figure left. 
The zincblende structure is displayed in figure^ middle, 
and has a similar energy as the wurtzite structure. At 
zero pressure, the energies of the wurtzite and zincblende 
structure are comparable to that of the layered structures 
in figure [T] 

The structure with space group 136 has six-membered 
rings (3 B, 3 N, alternating), but also rings with 4 (2 B, 2 
N) atoms. This leads to angles close to 90° and a less fa- 
vorable energy (see figure [2| right). This structure corre- 
sponds to the /3-BeO type (49|, which demonstrates that 
the method presented gives reasonable low-lying struc- 
ture candidates: as Be has one electron less than B, and 
O one more than N, it makes sense that such the /?- 
BeO structures is found as a candidate structure also for 
BN. Under ambient conditions, BeO crystallizes in the 
wurtzite structure, and the /3-BeO structure is found as 
a high temperature phase [49j. 

The structure with space group 62 has four-, six- and 
eight membered rings. The topology is similar to the one 
of the aluminum network in the SrAl 2 structure under 
ambient pressure (space group 74, Imma): one has to 
replace one aluminum atom with boron, the neighboring 
one with nitrogen, and discard the strontium. This is 
reasonable, as the two aluminum atoms obtain two elec- 
trons from strontium and thus have 8 valence electrons 
together, i.e. the same number of valence electrons as 
one boron and one nitrogen atom together. The struc- 
ture with space group 14 consists of layers - each of them 
consisting of rings with four or eight atoms. Finally, two 
structures with relatively large channels (i.e. large re- 
gions in the unit cell without atoms) were found, with 
space groups 8 and 9, respectively. 

The geometries are in reasonable agreement with the 
available experimental data in table HI The computed 
cell parameter a and the interlayer distance are approx- 
imatively constant for the layered structures with only 
six-membered rings (with space groups 194, 160, 187); 
this is also observed in the experiment when comparing 
the Bk and the rhombohedral structure. 

Total energies and statistics are given in table HVl The 



five energetically lowest lying structures found were the 
layered structures (space group 194, 160, 187) and the 
wurtzite and zincblende structure. At zero pressure, 
LDA favors wurtzite and zincblende (by 0.05 Eh per 
four formula units), whereas B3LYP (by 0.01 Eh) and 
Hartree-Fock favor the layered structures. This is in rea- 
sonable agreement with other calculations, e.g. reference 
[36| and references therein gives the layered structure by 
^0.06 eV/atom (0.02 Eh per four formula units) higher 
than the wurtzite and zincblende structures. Also, the 
zincblende and wurtzite structure are nearly degenerate, 
and similarly the hexagonal and rhombohedral structure. 
These results are stable with respect to the choice of the 
basis set, as calculations with larger basis sets show (see 
appendix [AJ . 

The statistics shows that the layered structures are 
frequently found, as well as the 3d-structures such as 
wurtzite or zincblende. The statistics includes all runs 
where Hartree-Fock energies were used during the sim- 
ulated annealing. 10 runs were performed, where the 
B3LYP functional was used during the simulated anneal- 
ing procedure. In two of these runs, a good structure 
candidate was found (wurtzite and the layered structure 
with space group 160). However, as was mentioned, the 
B3LYP runs require much more CPU time, and there- 
fore the runs were mainly performed using the Hartree- 
Fock approach during the simulated annealing. In total, 
about 16 % of the runs gave one of the structure can- 
didates in table IIV1 The other runs yielded either no 
good structure candidates (like " amorphous" structures) 
or only energetically very unfavorable structures. 

When pressure is applied, the lower coordinated struc- 
tures become less favorable, which is in agreement with 
the rule that the coordination number increases with 
pressure, see e.g. [13 ■ The enthalpies are displayed in 
figure [5] for LDA, and for B3LYP, respectively. Interest- 
ingly, in the case of B3LYP the layered structures are 
favorable at zero pressure, with wurtzite and zincblende 
becoming favorable at a pressure of ~ 3 GPa. 



5 




FIG. 2: (Color online) The structures found, with space group 186, 216 and 136. 





FIG. 3: (Color online) The structures found, with space group 62 and 8. 





FIG. 5: (Color online) The enthalpy of the most relevant structures, at the LDA and B3LYP level. Structure candidates are 
labeled by their space group. 




TABLE III: Basis sets used for the global search (I, II) and the local optimization (III) 



basis set I 
exponent contraction 



basis set II 
exponent contraction 



basis set III 
exponent contraction 



2.082E+03 
3.123E+02 
7.089E+01 
1.985E+01 
6.292E+00 
2.129E+00 



1.850E-03 
1.413E-02 
6.927E-02 
2.324E-01 
4.702E-01 
3.603E-01 



2.282E+00 -3.687E-01 2.312E-01 



B 

s 



2.082E+03 
3.123E+02 
7.089E+01 
1.985E+01 
6.292E+00 
2.129E+00 



1.850E-03 
1.413E-02 
6.927E-02 
2.324E-01 
4.702E-01 
3.603E-01 



sp 



0.4 



1.0 1.0 



2.282E+00 -3.687E-01 2.312E-01 
4.652E-01 1.199E+00 8.668E-01 
sp 

0.25 1.0 1.0 



2.082E+03 
3.123E+02 
7.089E+01 
1.985E+01 
6.292E+00 
2.129E+00 



1.850E-03 
1.413E-02 
6.927E-02 
2.324E-01 
4.702E-01 
3.603E-01 



2.282E+00 -3.687E-01 2.312E-01 
4.652E-01 1.199E+00 8.668E-01 



0.197 



1.0 1.0 



4150.0 
620.1 
141.7 
40.34 
13.03 
4.47 

5.425 
1.149 

0.3 



0.001845 
0.01416 
0.06863 
0.2286 
0.4662 
0.3657 

-0.4133 0.238 
1.224 0.859 

1.0 1.0 



4150.0 
620.1 
141.7 
40.34 
13.03 
4.47 

5.425 
1.149 

0.297 



N 

s 

0.001845 
0.01416 
0.06863 
0.2286 
0.4662 
0.3657 

sp 

-0.4133 0.238 
1.224 0.859 

sp 

1.0 1.0 

d 



0.8 

4150.0 
620.1 
141.7 
40.34 
13.03 
4.47 

5.425 
1.149 

0.297 

0.8 



1.0 



0.001845 
0.01416 
0.06863 
0.2286 
0.4662 
0.3657 

-0.4133 0.238 
1.224 0.859 

1.0 1.0 

1.0 



8 



TABLE IV: Total energies of the most relevant structures, and statistics. Energies are in hartree units (1 £^=27.2114 eV), for 
4 formula units. A run is considered successful, if one of the most relevant structures, as displayed in this table, was found. 



name of modification 


space group 


energy [E h 


1 


number of times found 






LDA 


B3LYP 


HF 


basis I 


basis II 


hexagonal BN 


194 


-315.9121 


-318.5115 


-316.8753 


1 


2 


I-BN 


160 


-315.9125 


-318.5110 


-316.8740 


11 


4 


II-BN 


187 


-315.9123 


-318.5109 


-316.8739 





2 


wurtzite 


186 


-315.9629 


-318.4988 


-316.8452 


8 


2 


zincblende 


216 


-315.9619 


-318.4991 


-316.8459 


4 


2 


/3-BeO 


136 


-315.9347 


-318.4753 


-316.8147 


2 


1 


III-BN 


62 


-315.8958 


-318.4437 


-316.7821 


4 


2 


IV-BN 


8 


-315.8810 


-318.4685 


-316.8243 


6 





V-BN 


9 


-315.8707 


-318.4569 


-316.8138 


2 





VI-BN 


14 


-315.8075 


-318.4085 


-316.7594 





1 



number of successful runs 36 (18.9 %) 16 (11.5 %) 

number of runs in total 190 139 



TABLE V: The energetically most favorable structures found. 



space group and 


cell parameters and fractional coordinates 


modification 


Ida 


RQT VP 


nr 


194 


a=2.50 A, c=5.88 A 


a=2.51 A, c=6.40 A 


a=2.50 A, c=6.43 A 


Bk 


B (1/3, 2/3, 1/4) 


B (1/3, 2/3, 1/4) 


B (1/3, 2/3, 1/4) 


hexagonal BN 


N (1/3, 2/3, 3/4) 


N (1/3, 2/3, 3/4) 


N (1/3, 2/3, 3/4) 


160 


a=2.50 A, c=8.72 A 


a= 2.51 A, c=9.56 A 


a=2.50 A, c=9.69 A 


I-BN 


D ^U, U, UJ 


a (u, u, uj 


a u, uj 




N (1/3, 2/3, -0.0009) 


N (1/3, 2/3, -0.0002) 


N (1/3, 2/3, 0.) 


187 


a=2.50 A, c=5.83 A 


a=2.51 A, c=6.38 A 


a=2.50 A, c=6.47 A 


II-BN 


a (U, u, uj 


a (U, U, UJ 


r fa n n\ 
a (U, U, UJ 




B (1/3, 2/3, 1/2) 


B (1/3, 2/3, 1/2) 


B (1/3, 2/3, 1/2) 




N (2/3, 1/3, 1/2) 


N (2/3, 1/3, 1/2) 


N (2/3, 1/3, 1/2) 




N (1/3, 2/3, 0) 


N (1/3, 2/3, 0) 


N (1/3, 2/3, 0) 


186 


a=2.54 A, c=4.19 A 


a=2.57 A, c=4.23 A 


a=2.55 A, c=4.21 A 


wurtzite 


B (2/3, 1/3, 0) 


B (2/3, 1/3, 0) 


B (2/3, 1/3, 0) 




N (2/3, 1/3, 0.3748) 


(2/3, 1/3, 0.3750) 


(2/3, 1/3, 0.3752) 


216 


a=3.60 A 


a=3.64 A 


a=3.62 A 


zincblende 


B (0, 0, 0) 


B (0, 0, 0) 


B (0, 0, 0) 




N (1/4, 1/4, 1/4) 


N (1/4, 1/4, 1/4) 


N (1/4, 1/4, 1/4) 


136 


a=4.38 A, c=2.54 A 


a=4.43 A, c=2.56 A 


a=4.41 A, c=2.55 A 


/3-BeO 


B (-0.1738, 0.1738, 1/2) 


B (-0.1744, 0.1744, 1/2) 


(-0.1742, 0.1742, 1/2) 




N (-0.1880, -0.1880, 1/2) 


N (-0.1872, -0.1872, 1/2) 


(-0.1872, -0.1872, 1/2) 


62 


a=4.76 A, b=2.58 A, c=4.29 A 


a=4.85 A, b=2.60 A, c=4.31 A 


a=4.83 A, b=2.59 A, c=4.30 A 


III-BN 


B (-0.3404, 3/4, 0.0921) 


B (-0.3397, 3/4, 0.0912) 


B (-0.3399, 3/4, 0.0905) 




N (0.3171, 3/4, 0.1048) 


N (0.3211, 3/4, 0.1074) 


N (0.3215, 3/4, 0.1078) 


8 


a=12.95 A, b=2.50 A, c=4.33 A 


a=13.08 A, b=2.51 A, c=4.39 A 


a=13.03 A, b=2.50 A, c=4.37 A 


IV-BN 


/?=91.7 


/3=90.7° 


,3=90.7° 




a (u, o, uj 


d ^U.UUUlj U, -U.UUllj 


d ^U.UUUz, U, -U.UUUoJ 




B (0.1998, 0., 0.4266) 


B (0.2000, 0., 0.4268) 


B (0.2003, 0., 0.4265) 




B (-0.1342, 0., 0.3784) 


B (-0.1344, 0., 0.3806) 


B (-0.1343, 0., 0.3820) 




B (0.0319, 1/2, -0.4995) 


B (0.0319, 1/2, 0.4995) 


B (0.0320, 1/2, 0.4998) 




N (-0.0192, 0., 0.3442) 


N (-0.0192, 0., 0.3436) 


N (-0.0194, 0., 0.3430) 




N (0.0062, 1/2, -0.1562) 


N (0.0057, 1/2, -0.1562) 


N (0.0051, 1/2, -0.1563) 




N (0.1468, 1/2, 0.4443) 


N (0.1471, 1/2, 0.4424) 


N (0.1471, 1/2, 0.4409) 




1ST ( fl 1 QTQ 1 /O f) QOflQ^ 

IN ^-U.iofO, 1/^, U.oyUoJ 


In ^-U.lo/o, 1/^, U.oyzoJ 


IN ^-U.lofO, 1/^, U.oyo4J 


9 


a=3.00 A, b=9.48 A, c=4.33 A 


a=3.21 A, b=9.46 A, c=4.35 A 


a=3.24 A, b=9.42 A, c=4.33 A 


V-BN 


(3=105.1° 


,3=105.8° 


P= 106.0° 




B (-0.4996, -0.0652, 0.0001) 


B (-0.4991, -0.0657, 0.0004) 


B (-0.4985, -0.0658, 0.0006) 




B (0.3780, -0.1931, 0.4920) 


B (0.3717, -0.1931, 0.4920) 


B (0.3711, -0.1929, 0.4923) 




N (-0.4219, -0.1948, -0.1601) 


N (-0.4173, -0.1955, -0.1593) 


N (-0.4160, -0.1953, -0.1592) 




N (0.4422, -0.0651, 0.3218) 


N (0.4435, -0.0654, 0.3209) 


N (0.4423, -0.0653, 0.3202) 


14 


a=3.17 A, b=4.90 A, c=4.93 


a=3.44 A, b=4.92 A, c=4.94 A, 


a=3.48 A, b=4.90 A, c=4.91 A 


VI-BN 


(3=117.7° 


(3=115.1° 


(3=114.3° 




B (0.4963, -0.1429, 0.1385) 


B (-0.4998, -0.1428, 0.1413) 


B (-0.4992, -0.1427, 0.1417) 




N (0.4984, 0.3443 -0.3418) 


N (0.4996, 0.3437, -0.3427) 


N (0.4989, 0.3441, -0.3437) 
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IV. CONCLUSION 

It was shown that structure prediction based on simu- 
lated annealing and using ab-initio energies during both 
the global and local optimization is feasible for a cova- 
lent system such as boron nitride. This is a significant 
extension of the previous work where this approach 
was shown to be feasible for an ionic system. Covalent 
systems are more difficult to study as covalent bonds 
need to be established between the neighbors, and con- 
vergence problems are more severe in this case. Three 
layered structures, the wurtzite and zincblende structure, 
a structure of the /3-BeO type, and four other favorable 
structures were found. Applying pressure leads to a pref- 
erence of the higher coordinated structures. 

APPENDIX A: LARGER BASIS SETS TO 
EXTRAPOLATE THE BASIS SET LIMIT 

In the present work, a main task is to compute the en- 
ergy differences between various structures. Besides the 
functional, the choice of the basis set has an influence on 
these results. In order to investigate this in more detail, 
the basis set used for the local optimization was further 
enlarged, and the total energies were computed for the 



most important structures. Here, enlarging the basis set 
means to include more diffuse functions. It turned out 
that this was only possible for the nitrogen atom, whereas 
more diffuse functions on boron led to linear dependence 
problems. Therefore, in a first step, one sp shell with 
exponent 0.15 was added to the nitrogen basis set III in 
table Hill which resulted in a [4s3pld] basis set for nitro- 
gen (basis set IV). In a second step, two sp shells (with 
exponents 0.15 and 0.6) were added to the nitrogen basis, 
i.e. a [5s4pld] basis set was obtained (basis set V). Note 
that these basis sets work reasonably well at zero pres- 
sure, but numerical instability sets in with compression, 
i.e. enthalpies (as in figure [5]) could only be obtained up 
to a relatively small pressure. 

A geometry optimization was performed with these ba- 
sis sets. The results are displayed in table I VII It becomes 
obvious, that the energy differences between the various 
structures remain essentially constant; the total energy 
becomes lower with increasing basis set. This is visual- 
ized in figure El 

The geometry slightly changes when enlarging the ba- 
sis set. Most prominent is the change of the c-axis for the 
layered structures when the basis set is enlarged. This is 
due to the weak bonding between the individual layers. 
However, as a whole, the basis set does not change the 
relative energies between the structures. 
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TABLE VI: A comparison of three basis sets for the energet- 
ically most favorable structures found. 



space group basis set cell parameters, in A energy, in E h 
and modification (LDA) 



194 


III 


a=2.50 c=5.88 


-315.9121 


hexagonal BN 


IV 


a=2.51 c=6.24 


-315.9556 




V 


a=2.51 c=6.18 


-315.9659 


160 


III 


a=2.50 c=8.72 


-315.9125 


I-BN 


IV 


a=2.51 c=9.27 


-315.9562 




V 


a=2.50 c=9.19 


-315.9668 


187 


III 


a=2.50 c=5.83 


-315.9123 


II-BN 


IV 


a=2.51 c=6.18 


-315.9564 




V 


a=2.50 c=6.12 


-315.9672 


186 


III 


a=2.54 c=4.19 


-315.9629 


wurtzite 


IV 


a=2.54 c=4.19 


-315.9884 




V 


a=2.54 c=4.18 


-316.0035 


ZIO 


TTT 
111 


a — o.du 


-oio.yoiy 


zincblende 


IV 


a=3.60 


-315.9873 




V 


a=3.60 


-316.0014 


136 


III 


a=4.38 c=2.54 


-315.9347 


III-BN 


IV 


a=4.39 c=2.54 


-315.9591 




V 


a=4.37 c=2.54 


-315.9749 


62 


III 


a=4.76 b=2.58 c=4.29 


-315.8958 


IV-BN 


IV 


a=4.79 b=2.58 c=4.29 


-315.9216 




V 


a=4.76 b=2.57 c=4.29 


-315.9366 
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FIG. 6: Total energies for the various structures, in hartree 
per four formula units, for basis sets III, IV and V, at the 
geometry optimized for each basis set. Structure candidates 
are labeled by the space group 



-, 1 1 1 1 1 1 — 

X basis set III 

+ basis set III and N sp 0.15 

O basis set III and Nsp 0.6, 0.15 



± + + V X + 

O O O xx 



+ + 

o o 



o 



+ 

o 



187 186 216 

structures (space group) 
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